/*


!
!  Dalton, a molecular electronic structure program
!  Copyright (C) The Dalton Authors (see AUTHORS file for details).
!
!  This program is free software; you can redistribute it and/or
!  modify it under the terms of the GNU Lesser General Public
!  License version 2.1 as published by the Free Software Foundation.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  Lesser General Public License for more details.
!
!  If a copy of the GNU LGPL v2.1 was not distributed with this
!  code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
!

!

*/
/*-*-mode: C; c-indentation-style: "bsd"; c-basic-offset: 4; -*-*/
/* fun-ft97bx.c:

   Automatically generated code implementing FT97BX functional and
   its derivatives. It is generated by func-codegen.pl being a part of
   a "Automatic code generation framework for analytical functional
   derivative evaluation", Pawel Salek, 2005

    This functional is connected by making following changes:
    1. add "extern Functional ft97bxFunctional;" to 'functionals.h'
    2. add "&ft97bxFunctional," to 'functionals.c'
    3. add "fun-ft97bx.c" to 'Makefile.am', 'Makefile.in' or 'Makefile'.

    This functional has been generated from following input:
    ------ cut here -------
rho:  rhoa + rhob;
grad: sqrt(grada*grada + gradb*gradb + 2*gradab);
zeta: (rhoa-rhob)/(rhoa+rhob);

bet0: 0.002913644;
bet1: 0.000947417;
bet2: 2501.149;
beta:  bet0+bet1*grada^2/(bet2+grada^2);
betb:  bet0+bet1*gradb^2/(bet2+gradb^2);

xa:   grada*rhoa^(-4/3);
xb:   gradb*rhob^(-4/3);

sa:   asinh(xa^2);
sb:   asinh(xb^2);

da:   1+9*beta^2*xa^2*sa^2;
db:   1+9*betb^2*xb^2*sb^2;

Exa:  -rhoa^(4/3)*beta*xa^2/(da^(1/2));
Exb:  -rhob^(4/3)*betb*xb^2/(db^(1/2));

K(rhoa,rhob,grada,gradb,gradab):=(Exa+Exb);


    ------ cut here -------
*/

 
/* strictly conform to XOPEN ANSI C standard */
#if !defined(SYS_DEC)
/* XOPEN compliance is missing on old Tru64 4.0E Alphas and pow() prototype
 * is not specified. */
#define _XOPEN_SOURCE          500
#define _XOPEN_SOURCE_EXTENDED 1
#endif
#include <math.h>
#include <stddef.h>
#include "general.h"

#define __CVERSION__
 
#include "functionals.h"
 
/* INTERFACE PART */
static integer ft97bx_isgga(void) { return 1; } /* FIXME: detect! */
static integer ft97bx_read(const char *conf_line);
static real ft97bx_energy(const FunDensProp* dp);
static void ft97bx_first(FunFirstFuncDrv *ds,   real factor,
                         const FunDensProp* dp);
static void ft97bx_second(FunSecondFuncDrv *ds, real factor,
                          const FunDensProp* dp);
static void ft97bx_third(FunThirdFuncDrv *ds,   real factor,
                         const FunDensProp* dp);
static void ft97bx_fourth(FunFourthFuncDrv *ds,   real factor,
                          const FunDensProp* dp);
 
Functional FT97bxFunctional = {
  "FT97bx",       /* name */
  ft97bx_isgga,   /* gga-corrected */
   1,
  ft97bx_read,
  NULL,
  ft97bx_energy,
  ft97bx_first,
  ft97bx_second,
  ft97bx_third,
  ft97bx_fourth
};
 
/* IMPLEMENTATION PART */
static integer
ft97bx_read(const char *conf_line)
{
    fun_set_hf_weight(0);
    return 1;
}

static real
ft97bx_energy(const FunDensProp *dp)
{
    real res;
    real rhoa = dp->rhoa, rhob = dp->rhob;
    real grada = dp->grada, gradb = dp->gradb, gradab = dp->gradab;

    real t1, t2, t3, t4, t5, t6;

    t1 = pow(grada,2.0);
    t2 = 9.47417E-4*t1/(t1+2501.149)+0.002913644;
    t3 = 1/pow(rhoa,2.666666666666667);
    t4 = pow(gradb,2.0);
    t5 = 9.47417E-4*t4/(t4+2501.149)+0.002913644;
    t6 = 1/pow(rhob,2.666666666666667);

   /* code */
    res = -1.0*t4*t5/(sqrt(9.0*t4*pow(t5,2.0)*t6*pow(asinh(t4*
        t6),2.0)+1.0)*pow(rhob,1.333333333333333))-1.0*t1*t2/(sqrt(9.0*
        t1*pow(t2,2.0)*t3*pow(asinh(t1*t3),2.0)+1.0)*pow(rhoa,1.333333333333333));

    return res;
}

static void
ft97bx_first_helper(real rhoa, real grada, real *res)
{    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17;

    t1 = pow(grada,2.0);
    t2 = t1+2501.149;
    t3 = 1/t2;
    t4 = 9.47417E-4*t1*t3+0.002913644;
    t5 = pow(t4,2.0);
    t6 = 1/pow(rhoa,2.666666666666667);
    t7 = asinh(t1*t6);
    t8 = pow(t7,2.0);
    t9 = sqrt(9.0*t1*t5*t6*t8+1.0);
    t10 = 1/t9;
    t11 = pow(grada,4.0);
    t12 = 1/pow(rhoa,5.333333333333333);
    t13 = 1/sqrt(t11*t12+1.0);
    t14 = 1/pow(t9,3.0);
    t15 = 1/pow(rhoa,1.333333333333333);
    t16 = pow(grada,3.0);
    t17 = 0.001894834*grada*t3-0.001894834*t16/pow(t2,2.0);

   /*
         code */
    res[0] = 0.5*t1*t14*t15*t4*(-48.0*t11*t13*t5*t7/pow(rhoa,
        6.333333333333333)-24.0*t1*t5*t8/pow(rhoa,3.666666666666667))+
        1.333333333333333*t1*t10*t4/pow(rhoa,2.333333333333333);
    res[1] = 0.5*t1*t14*t15*t4*(18.0*t5*t6*t8*grada+18.0*
        t1*t17*t4*t6*t8+36.0*t12*t13*t16*t5*t7)-2.0*t10*t15*t4*grada-
        1.0*t1*t10*t15*t17;
}

static void
ft97bx_first(FunFirstFuncDrv *ds, real factor, const FunDensProp *dp)
{
    real res[2];

    ft97bx_first_helper(dp->rhoa, dp->grada, res);
   /* Final assignment */
    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        ft97bx_first_helper(dp->rhob, dp->gradb, res);
    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

}

static void
ft97bx_second_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20, t21, t22, t23, t24, t25, t26;
    real t27, t28, t29, t30;

    t1 = pow(grada,2.0);
    t2 = t1+2501.149;
    t3 = 1/t2;
    t4 = 9.47417E-4*t1*t3+0.002913644;
    t5 = pow(t4,2.0);
    t6 = 1/pow(rhoa,2.666666666666667);
    t7 = asinh(t1*t6);
    t8 = pow(t7,2.0);
    t9 = sqrt(9.0*t1*t5*t6*t8+1.0);
    t10 = 1/t9;
    t11 = 1/pow(rhoa,2.333333333333333);
    t12 = pow(grada,4.0);
    t13 = 1/pow(rhoa,5.333333333333333);
    t14 = t12*t13+1.0;
    t15 = sqrt(t14);
    t16 = 1/t15;
    t17 = 1/pow(rhoa,6.333333333333333);
    t18 = 1/pow(rhoa,3.666666666666667);
    t19 = -24.0*t1*t18*t5*t8-48.0*t12*t16*t17*t5*t7;
    t20 = 1/pow(t9,3.0);
    t21 = 1/pow(rhoa,1.333333333333333);
    t22 = pow(grada,3.0);
    t23 = 1/pow(t2,2.0);
    t24 = 0.001894834*grada*t3-0.001894834*t22*t23;
    t25 = 18.0*t5*t6*t8*grada+18.0*t1*t24*t4*t6*t8+36.0*t13*
        t16*t22*t5*t7;
    t26 = 1/pow(t9,5.0);
    t27 = 1/pow(t15,3.0);
    t28 = pow(grada,6.0);
    t29 = 1/t14;
    t30 = 0.001894834*t3-0.00947417*t1*t23+0.007579336*t12/
        pow(t2,3.0);

   /* code */
    res[0] = 0.5*t1*t19*t20*t21*t4+1.333333333333333*t1*t10*
        t11*t4;
    res[1] = -2.0*t10*t21*t4*grada+0.5*t1*t20*t21*t25*t4-
        1.0*t1*t10*t21*t24;
    res[2] = 0.5*t1*t20*t21*t4*(-128.0*t27*t5*t7*pow(grada,
        8.0)/pow(rhoa,12.66666666666667)+128.0*t28*t29*t5/pow(rhoa,
        10.0)+432.0*t12*t16*t5*t7/pow(rhoa,7.333333333333333)+88.0*
        t1*t5*t8/pow(rhoa,4.666666666666667))-3.111111111111111*t1*
        t10*t4/pow(rhoa,3.333333333333333)-0.75*t1*pow(t19,2.0)*t21*
        t26*t4-1.333333333333333*t1*t11*t19*t20*t4;
    res[3] = 0.5*t1*t20*t21*t4*(96.0*t27*t5*t7*pow(grada,
        7.0)/pow(rhoa,11.66666666666667)-96.0*t29*t5*pow(grada,5.0)/
        pow(rhoa,9.0)-48.0*t18*t5*t8*grada-48.0*t1*t18*t24*t4*t8-288.0*
        t16*t17*t22*t5*t7-96.0*t12*t16*t17*t24*t4*t7)+2.666666666666667*
        t10*t11*t4*grada-0.75*t1*t19*t21*t25*t26*t4-0.666666666666667*
        t1*t11*t20*t25*t4+0.5*t1*t19*t20*t21*t24+1.333333333333333*
        t1*t10*t11*t24+grada*t4*t19*t20*t21;
    res[4] = 0.5*t1*t20*t21*t4*(-72.0*t27*t28*t5*t7/pow(rhoa,
        10.66666666666667)+72.0*t12*t29*t5/pow(rhoa,8.0)+72.0*t24*
        t4*t6*t8*grada+18.0*t5*t6*t8+18.0*t1*t30*t4*t6*t8+18.0*t1*
        pow(t24,2.0)*t6*t8+180.0*t1*t13*t16*t5*t7+144.0*t13*t16*t22*
        t24*t4*t7)+2.0*t20*t21*t25*t4*grada-4.0*t10*t21*t24*grada-
        0.75*t1*t21*pow(t25,2.0)*t26*t4-2.0*t10*t21*t4-1.0*t1*t10*
        t21*t30+t1*t24*t25*t20*t21;

}

static void
ft97bx_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[5];
 
    ft97bx_second_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        ft97bx_second_helper(dp->rhob, dp->gradb, res);
    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

}

static void
ft97bx_third_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20, t21, t22, t23, t24, t25, t26;
    real t27, t28, t29, t30, t31, t32, t33, t34;
    real t35, t36, t37, t38, t39, t40, t41, t42;
    real t43, t44, t45, t46, t47, t48, t49, t50;
    real t51, t52, t53, t54, t55;

    t1 = pow(grada,2.0);
    t2 = t1+2501.149;
    t3 = 1/t2;
    t4 = 9.47417E-4*t1*t3+0.002913644;
    t5 = pow(t4,2.0);
    t6 = 1/pow(rhoa,2.666666666666667);
    t7 = asinh(t1*t6);
    t8 = pow(t7,2.0);
    t9 = sqrt(9.0*t1*t5*t6*t8+1.0);
    t10 = 1/t9;
    t11 = 1/pow(rhoa,2.333333333333333);
    t12 = pow(grada,4.0);
    t13 = 1/pow(rhoa,5.333333333333333);
    t14 = t12*t13+1.0;
    t15 = sqrt(t14);
    t16 = 1/t15;
    t17 = 1/pow(rhoa,6.333333333333333);
    t18 = 1/pow(rhoa,3.666666666666667);
    t19 = -24.0*t1*t18*t5*t8-48.0*t12*t16*t17*t5*t7;
    t20 = 1/pow(t9,3.0);
    t21 = 1/pow(rhoa,1.333333333333333);
    t22 = pow(grada,3.0);
    t23 = 1/pow(t2,2.0);
    t24 = 0.001894834*grada*t3-0.001894834*t22*t23;
    t25 = 18.0*t5*t6*t8*grada+18.0*t1*t24*t4*t6*t8+36.0*t13*
        t16*t22*t5*t7;
    t26 = 1/pow(rhoa,3.333333333333333);
    t27 = pow(t19,2.0);
    t28 = 1/pow(t9,5.0);
    t29 = pow(grada,8.0);
    t30 = 1/pow(t15,3.0);
    t31 = 1/pow(rhoa,12.66666666666667);
    t32 = pow(grada,6.0);
    t33 = 1/t14;
    t34 = 1/pow(rhoa,10.0);
    t35 = 1/pow(rhoa,7.333333333333333);
    t36 = 1/pow(rhoa,4.666666666666667);
    t37 = 88.0*t1*t36*t5*t8+432.0*t12*t16*t35*t5*t7-128.0*
        t29*t30*t31*t5*t7+128.0*t32*t33*t34*t5;
    t38 = pow(grada,7.0);
    t39 = 1/pow(rhoa,11.66666666666667);
    t40 = pow(grada,5.0);
    t41 = 1/pow(rhoa,9.0);
    t42 = -48.0*t18*t5*t8*grada-48.0*t1*t18*t24*t4*t8+96.0*
        t30*t38*t39*t5*t7-288.0*t16*t17*t22*t5*t7-96.0*t12*t16*t17*
        t24*t4*t7-96.0*t33*t40*t41*t5;
    t43 = pow(t25,2.0);
    t44 = 1/pow(rhoa,10.66666666666667);
    t45 = 1/pow(rhoa,8.0);
    t46 = pow(t24,2.0);
    t47 = 1/pow(t2,3.0);
    t48 = 0.001894834*t3-0.00947417*t1*t23+0.007579336*t12*
        t47;
    t49 = 72.0*t24*t4*t6*t8*grada+18.0*t5*t6*t8+18.0*t1*t4*
        t48*t6*t8+18.0*t1*t46*t6*t8-72.0*t30*t32*t44*t5*t7+180.0*t1*
        t13*t16*t5*t7+144.0*t13*t16*t22*t24*t4*t7+72.0*t12*t33*t45*
        t5;
    t50 = 1/pow(t9,7.0);
    t51 = 1/pow(t15,5.0);
    t52 = pow(grada,10.0);
    t53 = 1/pow(t14,2.0);
    t54 = pow(grada,9.0);
    t55 = 0.068214024*t22*t47-0.045476016*t40/pow(t2,4.0)-
        0.022738008*grada*t23;

   /* code */
    res[0] = 0.5*t1*t19*t20*t21*t4+1.333333333333333*t1*t10*
        t11*t4;
    res[1] = -2.0*t10*t21*t4*grada+0.5*t1*t20*t21*t25*t4-
        1.0*t1*t10*t21*t24;
    res[2] = 0.5*t1*t20*t21*t37*t4-0.75*t1*t21*t27*t28*t4-
        3.111111111111111*t1*t10*t26*t4-1.333333333333333*t1*t11*t19*
        t20*t4;
    res[3] = 2.666666666666667*t10*t11*t4*grada+0.5*t1*t20*
        t21*t4*t42-0.75*t1*t19*t21*t25*t28*t4-0.666666666666667*t1*
        t11*t20*t25*t4+0.5*t1*t19*t20*t21*t24+1.333333333333333*t1*
        t10*t11*t24+grada*t4*t19*t20*t21;
    res[4] = 2.0*t20*t21*t25*t4*grada-4.0*t10*t21*t24*grada+
        0.5*t1*t20*t21*t4*t49-1.0*t1*t10*t21*t48-0.75*t1*t21*t28*t4*
        t43-2.0*t10*t21*t4+t1*t24*t25*t20*t21;
    res[5] = 0.5*t1*t20*t21*t4*(-1024.0*t5*t51*t7*pow(grada,
        12.0)/pow(rhoa,19.0)+1024.0*t5*t52*t53/pow(rhoa,16.33333333333333)+
        2773.333333333333*t29*t30*t5*t7/pow(rhoa,13.66666666666667)-
        2432.0*t32*t33*t5/pow(rhoa,11.0)-3637.333333333333*t12*t16*
        t5*t7/pow(rhoa,8.333333333333334)-410.6666666666667*t1*t5*
        t8/pow(rhoa,5.666666666666667))+10.37037037037037*t1*t10*t4/
        pow(rhoa,4.333333333333333)+1.875*t1*pow(t19,3.0)*t21*t4*t50-
        2.25*t1*t19*t21*t28*t37*t4-2.0*t1*t11*t20*t37*t4+3.0*t1*t11*
        t27*t28*t4+4.666666666666667*t1*t19*t20*t26*t4;
    res[6] = 0.5*t1*t20*t21*t4*(768.0*t5*t51*t7*pow(grada,
        11.0)/pow(rhoa,18.0)-768.0*t5*t53*t54/pow(rhoa,15.33333333333333)+
        176.0*t36*t5*t8*grada+176.0*t1*t24*t36*t4*t8-1888.0*t30*t31*
        t38*t5*t7+2080.0*t16*t22*t35*t5*t7+864.0*t12*t16*t24*t35*t4*
        t7-256.0*t24*t29*t30*t31*t4*t7+1632.0*t33*t34*t40*t5+256.0*
        t24*t32*t33*t34*t4)-1.5*t21*t27*t28*t4*grada-6.222222222222222*
        t10*t26*t4*grada-2.666666666666667*t11*t19*t20*t4*grada+1.875*
        t1*t21*t25*t27*t4*t50-1.5*t1*t19*t21*t28*t4*t42-1.333333333333333*
        t1*t11*t20*t4*t42-0.75*t1*t21*t25*t28*t37*t4+2.0*t1*t11*t19*
        t25*t28*t4+1.555555555555556*t1*t20*t25*t26*t4+0.5*t1*t20*
        t21*t24*t37-0.75*t1*t21*t24*t27*t28-3.111111111111111*t1*t10*
        t24*t26-1.333333333333333*t1*t11*t19*t20*t24+grada*t4*t37*
        t20*t21;
    res[7] = 0.5*t1*t20*t21*t4*(-576.0*t5*t51*t52*t7/pow(rhoa,
        17.0)+576.0*t29*t5*t53/pow(rhoa,14.33333333333333)-192.0*t18*
        t24*t4*t8*grada-48.0*t18*t5*t8-48.0*t1*t18*t4*t48*t8-48.0*
        t1*t18*t46*t8+1248.0*t30*t32*t39*t5*t7-1056.0*t1*t16*t17*t5*
        t7-96.0*t12*t16*t17*t4*t48*t7-96.0*t12*t16*t17*t46*t7+384.0*
        t24*t30*t38*t39*t4*t7-1152.0*t16*t17*t22*t24*t4*t7-1056.0*
        t12*t33*t41*t5-384.0*t24*t33*t4*t40*t41)+2.0*t20*t21*t4*t42*
        grada-3.0*t19*t21*t25*t28*t4*grada-2.666666666666667*t11*t20*
        t25*t4*grada+2.0*t19*t20*t21*t24*grada+5.333333333333333*t10*
        t11*t24*grada+1.875*t1*t19*t21*t4*t43*t50-0.75*t1*t19*t21*
        t28*t4*t49-0.666666666666667*t1*t11*t20*t4*t49+0.5*t1*t19*
        t20*t21*t48+1.333333333333333*t1*t10*t11*t48-1.5*t1*t21*t25*
        t28*t4*t42+2.666666666666667*t10*t11*t4-1.5*t1*t19*t21*t24*
        t25*t28-1.333333333333333*t1*t11*t20*t24*t25+t1*t24*t42*t20*
        t21+t4*t19*t20*t21+t1*t4*t43*t28*t11;
    res[8] = 0.5*t1*t20*t21*t4*(432.0*t5*t51*t54*t7/pow(rhoa,
        16.0)-432.0*t38*t5*t53/pow(rhoa,13.33333333333333)+108.0*t4*
        t48*t6*t8*grada+108.0*t46*t6*t8*grada+432.0*t13*t16*t5*t7*
        grada+18.0*t1*t4*t55*t6*t8+54.0*t1*t24*t48*t6*t8+108.0*t24*
        t4*t6*t8-792.0*t30*t40*t44*t5*t7+216.0*t13*t16*t22*t4*t48*
        t7+216.0*t13*t16*t22*t46*t7-432.0*t24*t30*t32*t4*t44*t7+1080.0*
        t1*t13*t16*t24*t4*t7+648.0*t22*t33*t45*t5+432.0*t12*t24*t33*
        t4*t45)+3.0*t20*t21*t4*t49*grada-6.0*t10*t21*t48*grada-4.5*
        t21*t28*t4*t43*grada+6.0*t20*t21*t24*t25*grada-1.0*t1*t10*
        t21*t55+1.875*t1*t21*pow(t25,3.0)*t4*t50-2.25*t1*t21*t25*t28*
        t4*t49+1.5*t1*t20*t21*t24*t49+1.5*t1*t20*t21*t25*t48-2.25*
        t1*t21*t24*t28*t43+3.0*t20*t21*t25*t4-6.0*t10*t21*t24;

}

static void
ft97bx_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[9];
 
    ft97bx_third_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];

    ds->df3000 += factor*res[5];
    ds->df2010 += factor*res[6];
    ds->df1020 += factor*res[7];
    ds->df0030 += factor*res[8];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        ft97bx_third_helper(dp->rhob, dp->gradb, res);

    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

    ds->df0300 += factor*res[5];
    ds->df0201 += factor*res[6];
    ds->df0102 += factor*res[7];
    ds->df0003 += factor*res[8];

}

static void
ft97bx_fourth_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20, t21, t22, t23, t24, t25, t26;
    real t27, t28, t29, t30, t31, t32, t33, t34;
    real t35, t36, t37, t38, t39, t40, t41, t42;
    real t43, t44, t45, t46, t47, t48, t49, t50;
    real t51, t52, t53, t54, t55, t56, t57, t58;
    real t59, t60, t61, t62, t63, t64, t65, t66;
    real t67, t68, t69, t70, t71, t72, t73, t74;
    real t75, t76, t77, t78, t79, t80, t81, t82;
    real t83;

    t1 = pow(grada,2.0);
    t2 = t1+2501.149;
    t3 = 1/t2;
    t4 = 9.47417E-4*t1*t3+0.002913644;
    t5 = pow(t4,2.0);
    t6 = 1/pow(rhoa,2.666666666666667);
    t7 = asinh(t1*t6);
    t8 = pow(t7,2.0);
    t9 = sqrt(9.0*t1*t5*t6*t8+1.0);
    t10 = 1/t9;
    t11 = 1/pow(rhoa,2.333333333333333);
    t12 = pow(grada,4.0);
    t13 = 1/pow(rhoa,5.333333333333333);
    t14 = t12*t13+1.0;
    t15 = sqrt(t14);
    t16 = 1/t15;
    t17 = 1/pow(rhoa,6.333333333333333);
    t18 = 1/pow(rhoa,3.666666666666667);
    t19 = -24.0*t1*t18*t5*t8-48.0*t12*t16*t17*t5*t7;
    t20 = 1/pow(t9,3.0);
    t21 = 1/pow(rhoa,1.333333333333333);
    t22 = pow(grada,3.0);
    t23 = 1/pow(t2,2.0);
    t24 = 0.001894834*grada*t3-0.001894834*t22*t23;
    t25 = 18.0*t5*t6*t8*grada+18.0*t1*t24*t4*t6*t8+36.0*t13*
        t16*t22*t5*t7;
    t26 = 1/pow(rhoa,3.333333333333333);
    t27 = pow(t19,2.0);
    t28 = 1/pow(t9,5.0);
    t29 = pow(grada,8.0);
    t30 = 1/pow(t15,3.0);
    t31 = 1/pow(rhoa,12.66666666666667);
    t32 = pow(grada,6.0);
    t33 = 1/t14;
    t34 = 1/pow(rhoa,10.0);
    t35 = 1/pow(rhoa,7.333333333333333);
    t36 = 1/pow(rhoa,4.666666666666667);
    t37 = 88.0*t1*t36*t5*t8+432.0*t12*t16*t35*t5*t7-128.0*
        t29*t30*t31*t5*t7+128.0*t32*t33*t34*t5;
    t38 = pow(grada,7.0);
    t39 = 1/pow(rhoa,11.66666666666667);
    t40 = pow(grada,5.0);
    t41 = 1/pow(rhoa,9.0);
    t42 = -48.0*t18*t5*t8*grada-48.0*t1*t18*t24*t4*t8+96.0*
        t30*t38*t39*t5*t7-288.0*t16*t17*t22*t5*t7-96.0*t12*t16*t17*
        t24*t4*t7-96.0*t33*t40*t41*t5;
    t43 = pow(t25,2.0);
    t44 = 1/pow(rhoa,10.66666666666667);
    t45 = 1/pow(rhoa,8.0);
    t46 = pow(t24,2.0);
    t47 = 1/pow(t2,3.0);
    t48 = 0.001894834*t3-0.00947417*t1*t23+0.007579336*t12*
        t47;
    t49 = 72.0*t24*t4*t6*t8*grada+18.0*t5*t6*t8+18.0*t1*t4*
        t48*t6*t8+18.0*t1*t46*t6*t8-72.0*t30*t32*t44*t5*t7+180.0*t1*
        t13*t16*t5*t7+144.0*t13*t16*t22*t24*t4*t7+72.0*t12*t33*t45*
        t5;
    t50 = 1/pow(rhoa,4.333333333333333);
    t51 = pow(t19,3.0);
    t52 = 1/pow(t9,7.0);
    t53 = pow(grada,12.0);
    t54 = 1/pow(t15,5.0);
    t55 = 1/pow(rhoa,19.0);
    t56 = pow(grada,10.0);
    t57 = 1/pow(t14,2.0);
    t58 = 1/pow(rhoa,16.33333333333333);
    t59 = 1/pow(rhoa,13.66666666666667);
    t60 = 1/pow(rhoa,11.0);
    t61 = 1/pow(rhoa,8.333333333333334);
    t62 = 1/pow(rhoa,5.666666666666667);
    t63 = -410.6666666666667*t1*t5*t62*t8-3637.333333333333*
        t12*t16*t5*t61*t7+2773.333333333333*t29*t30*t5*t59*t7-1024.0*
        t5*t53*t54*t55*t7-2432.0*t32*t33*t5*t60+1024.0*t5*t56*t57*
        t58;
    t64 = pow(grada,11.0);
    t65 = 1/pow(rhoa,18.0);
    t66 = pow(grada,9.0);
    t67 = 1/pow(rhoa,15.33333333333333);
    t68 = 176.0*t36*t5*t8*grada+176.0*t1*t24*t36*t4*t8+768.0*
        t5*t54*t64*t65*t7-1888.0*t30*t31*t38*t5*t7+2080.0*t16*t22*
        t35*t5*t7+864.0*t12*t16*t24*t35*t4*t7-256.0*t24*t29*t30*t31*
        t4*t7-768.0*t5*t57*t66*t67+1632.0*t33*t34*t40*t5+256.0*t24*
        t32*t33*t34*t4;
    t69 = 1/pow(rhoa,17.0);
    t70 = 1/pow(rhoa,14.33333333333333);
    t71 = -192.0*t18*t24*t4*t8*grada-48.0*t18*t5*t8-48.0*
        t1*t18*t4*t48*t8-48.0*t1*t18*t46*t8+576.0*t29*t5*t57*t70-576.0*
        t5*t54*t56*t69*t7+1248.0*t30*t32*t39*t5*t7-1056.0*t1*t16*t17*
        t5*t7-96.0*t12*t16*t17*t4*t48*t7-96.0*t12*t16*t17*t46*t7+384.0*
        t24*t30*t38*t39*t4*t7-1152.0*t16*t17*t22*t24*t4*t7-1056.0*
        t12*t33*t41*t5-384.0*t24*t33*t4*t40*t41;
    t72 = pow(t25,3.0);
    t73 = 1/pow(rhoa,16.0);
    t74 = 1/pow(rhoa,13.33333333333333);
    t75 = 1/pow(t2,4.0);
    t76 = -0.022738008*grada*t23+0.068214024*t22*t47-0.045476016*
        t40*t75;
    t77 = 108.0*t4*t48*t6*t8*grada+108.0*t46*t6*t8*grada+
        432.0*t13*t16*t5*t7*grada+18.0*t1*t4*t6*t76*t8+54.0*t1*t24*
        t48*t6*t8+108.0*t24*t4*t6*t8-432.0*t38*t5*t57*t74+432.0*t5*
        t54*t66*t7*t73-792.0*t30*t40*t44*t5*t7+216.0*t13*t16*t22*t4*
        t48*t7+216.0*t13*t16*t22*t46*t7-432.0*t24*t30*t32*t4*t44*t7+
        1080.0*t1*t13*t16*t24*t4*t7+648.0*t22*t33*t45*t5+432.0*t12*
        t24*t33*t4*t45;
    t78 = 1/pow(t9,9.0);
    t79 = 1/pow(t15,7.0);
    t80 = pow(grada,14.0);
    t81 = 1/pow(t14,3.0);
    t82 = pow(grada,13.0);
    t83 = -0.636664224*t12*t75+0.295594104*t1*t47+0.363808128*
        t32/pow(t2,5.0)-0.022738008*t23;

   /* code */
    res[0] = 0.5*t1*t19*t20*t21*t4+1.333333333333333*t1*t10*
        t11*t4;
    res[1] = -2.0*t10*t21*t4*grada+0.5*t1*t20*t21*t25*t4-
        1.0*t1*t10*t21*t24;
    res[2] = 0.5*t1*t20*t21*t37*t4-0.75*t1*t21*t27*t28*t4-
        3.111111111111111*t1*t10*t26*t4-1.333333333333333*t1*t11*t19*
        t20*t4;
    res[3] = 2.666666666666667*t10*t11*t4*grada+0.5*t1*t20*
        t21*t4*t42-0.75*t1*t19*t21*t25*t28*t4-0.666666666666667*t1*
        t11*t20*t25*t4+0.5*t1*t19*t20*t21*t24+1.333333333333333*t1*
        t10*t11*t24+grada*t4*t19*t20*t21;
    res[4] = 2.0*t20*t21*t25*t4*grada-4.0*t10*t21*t24*grada+
        0.5*t1*t20*t21*t4*t49-1.0*t1*t10*t21*t48-0.75*t1*t21*t28*t4*
        t43-2.0*t10*t21*t4+t1*t24*t25*t20*t21;
    res[5] = 0.5*t1*t20*t21*t4*t63+1.875*t1*t21*t4*t51*t52+
        10.37037037037037*t1*t10*t4*t50-2.25*t1*t19*t21*t28*t37*t4-
        2.0*t1*t11*t20*t37*t4+3.0*t1*t11*t27*t28*t4+4.666666666666667*
        t1*t19*t20*t26*t4;
    res[6] = -1.5*t21*t27*t28*t4*grada-6.222222222222222*
        t10*t26*t4*grada-2.666666666666667*t11*t19*t20*t4*grada+0.5*
        t1*t20*t21*t4*t68+1.875*t1*t21*t25*t27*t4*t52-1.5*t1*t19*t21*
        t28*t4*t42-1.333333333333333*t1*t11*t20*t4*t42-0.75*t1*t21*
        t25*t28*t37*t4+2.0*t1*t11*t19*t25*t28*t4+1.555555555555556*
        t1*t20*t25*t26*t4+0.5*t1*t20*t21*t24*t37-0.75*t1*t21*t24*t27*
        t28-3.111111111111111*t1*t10*t24*t26-1.333333333333333*t1*
        t11*t19*t20*t24+grada*t4*t37*t20*t21;
    res[7] = 2.0*t20*t21*t4*t42*grada-3.0*t19*t21*t25*t28*
        t4*grada-2.666666666666667*t11*t20*t25*t4*grada+2.0*t19*t20*
        t21*t24*grada+5.333333333333333*t10*t11*t24*grada+0.5*t1*t20*
        t21*t4*t71+1.875*t1*t19*t21*t4*t43*t52-0.75*t1*t19*t21*t28*
        t4*t49-0.666666666666667*t1*t11*t20*t4*t49+0.5*t1*t19*t20*
        t21*t48+1.333333333333333*t1*t10*t11*t48-1.5*t1*t21*t25*t28*
        t4*t42+2.666666666666667*t10*t11*t4-1.5*t1*t19*t21*t24*t25*
        t28-1.333333333333333*t1*t11*t20*t24*t25+t1*t24*t42*t20*t21+
        t4*t19*t20*t21+t1*t4*t43*t28*t11;
    res[8] = 3.0*t20*t21*t4*t49*grada-6.0*t10*t21*t48*grada-
        4.5*t21*t28*t4*t43*grada+6.0*t20*t21*t24*t25*grada+0.5*t1*
        t20*t21*t4*t77-1.0*t1*t10*t21*t76+1.875*t1*t21*t4*t52*t72-
        2.25*t1*t21*t25*t28*t4*t49+1.5*t1*t20*t21*t24*t49+1.5*t1*t20*
        t21*t25*t48-2.25*t1*t21*t24*t28*t43+3.0*t20*t21*t25*t4-6.0*
        t10*t21*t24;
    res[9] = 0.5*t1*t20*t21*t4*(-13653.33333333333*t5*t7*
        t79*pow(grada,16.0)/pow(rhoa,25.33333333333333)+13653.33333333333*
        t5*t80*t81/pow(rhoa,22.66666666666667)+41642.66666666666*t5*
        t53*t54*t7/pow(rhoa,20.0)-37091.55555555555*t5*t56*t57/pow(rhoa,
        17.33333333333333)-47601.77777777778*t29*t30*t5*t7/pow(rhoa,
        14.66666666666667)+36451.55555555555*t32*t33*t5/pow(rhoa,12.0)+
        32501.33333333333*t12*t16*t5*t7/pow(rhoa,9.333333333333334)+
        2327.111111111111*t1*t5*t8/pow(rhoa,6.666666666666667))-6.5625*
        t1*pow(t19,4.0)*t21*t4*t78-3.0*t1*t19*t21*t28*t4*t63-2.666666666666667*
        t1*t11*t20*t4*t63-10.0*t1*t11*t4*t51*t52+11.25*t1*t21*t27*
        t37*t4*t52-20.74074074074074*t1*t19*t20*t4*t50-2.25*t1*t21*
        t28*pow(t37,2.0)*t4+12.0*t1*t11*t19*t28*t37*t4+9.333333333333334*
        t1*t20*t26*t37*t4-14.0*t1*t26*t27*t28*t4-44.93827160493827*
        t1*t10*t13*t4;
    res[10] = 0.5*t1*t20*t21*t4*(10240.0*t5*t7*t79*pow(grada,
        15.0)/pow(rhoa,24.33333333333333)-10240.0*t5*t81*t82/pow(rhoa,
        21.66666666666667)-821.3333333333334*t5*t62*t8*grada-821.3333333333334*
        t1*t24*t4*t62*t8-28928.0*t5*t54*t55*t64*t7-16192.0*t16*t22*
        t5*t61*t7-7274.666666666667*t12*t16*t24*t4*t61*t7+29461.33333333333*
        t30*t38*t5*t59*t7+5546.666666666667*t24*t29*t30*t4*t59*t7-
        2048.0*t24*t4*t53*t54*t55*t7+25514.66666666667*t5*t57*t58*
        t66-21866.66666666667*t33*t40*t5*t60-4864.0*t24*t32*t33*t4*
        t60+2048.0*t24*t4*t56*t57*t58)+3.75*t21*t4*t51*t52*grada+20.74074074074074*
        t10*t4*t50*grada-4.5*t19*t21*t28*t37*t4*grada-4.0*t11*t20*
        t37*t4*grada+6.0*t11*t27*t28*t4*grada+9.333333333333334*t19*
        t20*t26*t4*grada-6.5625*t1*t21*t25*t4*t51*t78-2.25*t1*t19*
        t21*t28*t4*t68-2.0*t1*t11*t20*t4*t68-0.75*t1*t21*t25*t28*t4*
        t63+0.5*t1*t20*t21*t24*t63+1.875*t1*t21*t24*t51*t52+5.625*
        t1*t21*t27*t4*t42*t52+5.625*t1*t19*t21*t25*t37*t4*t52-7.5*
        t1*t11*t25*t27*t4*t52-5.185185185185185*t1*t20*t25*t4*t50+
        10.37037037037037*t1*t10*t24*t50-2.25*t1*t21*t28*t37*t4*t42+
        6.0*t1*t11*t19*t28*t4*t42+4.666666666666667*t1*t20*t26*t4*
        t42+3.0*t1*t11*t25*t28*t37*t4-7.0*t1*t19*t25*t26*t28*t4-2.25*
        t1*t19*t21*t24*t28*t37-2.0*t1*t11*t20*t24*t37+3.0*t1*t11*t24*
        t27*t28+4.666666666666667*t1*t19*t20*t24*t26+grada*t4*t63*
        t20*t21;
    res[11] = 0.5*t1*t20*t21*t4*(-7680.0*t5*t7*t79*t80/pow(rhoa,
        23.33333333333333)+7680.0*t5*t53*t81/pow(rhoa,20.66666666666667)+
        704.0*t24*t36*t4*t8*grada+176.0*t36*t5*t8+176.0*t1*t36*t4*
        t48*t8+176.0*t1*t36*t46*t8+3072.0*t24*t4*t54*t64*t65*t7+19776.0*
        t5*t54*t56*t65*t7+6944.0*t1*t16*t35*t5*t7-17376.0*t30*t31*
        t32*t5*t7+864.0*t12*t16*t35*t4*t48*t7-256.0*t29*t30*t31*t4*
        t48*t7+864.0*t12*t16*t35*t46*t7-256.0*t29*t30*t31*t46*t7-7552.0*
        t24*t30*t31*t38*t4*t7+8320.0*t16*t22*t24*t35*t4*t7-3072.0*
        t24*t4*t57*t66*t67-17216.0*t29*t5*t57*t67+12320.0*t12*t33*
        t34*t5+256.0*t32*t33*t34*t4*t48+256.0*t32*t33*t34*t46+6528.0*
        t24*t33*t34*t4*t40)+2.0*t20*t21*t4*t68*grada+7.5*t21*t25*t27*
        t4*t52*grada-6.0*t19*t21*t28*t4*t42*grada-5.333333333333333*
        t11*t20*t4*t42*grada-3.0*t21*t25*t28*t37*t4*grada+8.0*t11*
        t19*t25*t28*t4*grada+6.222222222222222*t20*t25*t26*t4*grada+
        2.0*t20*t21*t24*t37*grada-3.0*t21*t24*t27*t28*grada-12.44444444444444*
        t10*t24*t26*grada-5.333333333333333*t11*t19*t20*t24*grada-
        6.5625*t1*t21*t27*t4*t43*t78-1.5*t1*t19*t21*t28*t4*t71-1.333333333333333*
        t1*t11*t20*t4*t71-1.5*t1*t21*t25*t28*t4*t68+1.875*t1*t21*t27*
        t4*t49*t52+1.875*t1*t21*t37*t4*t43*t52-5.0*t1*t11*t19*t4*t43*
        t52+7.5*t1*t19*t21*t25*t4*t42*t52+3.75*t1*t21*t24*t25*t27*
        t52-0.75*t1*t21*t28*t37*t4*t49+2.0*t1*t11*t19*t28*t4*t49+1.555555555555556*
        t1*t20*t26*t4*t49+0.5*t1*t20*t21*t37*t48-0.75*t1*t21*t27*t28*
        t48-3.111111111111111*t1*t10*t26*t48-1.333333333333333*t1*
        t11*t19*t20*t48-2.333333333333333*t1*t26*t28*t4*t43-1.5*t1*
        t21*t28*t4*pow(t42,2.0)+4.0*t1*t11*t25*t28*t4*t42-3.0*t1*t19*
        t21*t24*t28*t42-2.666666666666667*t1*t11*t20*t24*t42-1.5*t21*
        t27*t28*t4-6.222222222222222*t10*t26*t4-2.666666666666667*
        t11*t19*t20*t4-1.5*t1*t21*t24*t25*t28*t37+4.0*t1*t11*t19*t24*
        t25*t28+3.111111111111111*t1*t20*t24*t25*t26+t1*t24*t68*t20*
        t21+t4*t37*t20*t21;
    res[12] = 0.5*t1*t20*t21*t4*(5760.0*t5*t7*t79*t82/pow(rhoa,
        22.33333333333333)-5760.0*t5*t64*t81/pow(rhoa,19.66666666666667)-
        288.0*t18*t4*t48*t8*grada-288.0*t18*t46*t8*grada-2304.0*t16*
        t17*t5*t7*grada-48.0*t1*t18*t4*t76*t8-144.0*t1*t18*t24*t48*
        t8-288.0*t18*t24*t4*t8-96.0*t12*t16*t17*t4*t7*t76+11328.0*
        t38*t5*t57*t70+3456.0*t24*t29*t4*t57*t70-13248.0*t5*t54*t66*
        t69*t7-3456.0*t24*t4*t54*t56*t69*t7+9600.0*t30*t39*t40*t5*
        t7+576.0*t30*t38*t39*t4*t48*t7-1728.0*t16*t17*t22*t4*t48*t7-
        288.0*t12*t16*t17*t24*t48*t7+576.0*t30*t38*t39*t46*t7-1728.0*
        t16*t17*t22*t46*t7+7488.0*t24*t30*t32*t39*t4*t7-6336.0*t1*
        t16*t17*t24*t4*t7-6336.0*t22*t33*t41*t5-576.0*t33*t4*t40*t41*
        t48-576.0*t33*t40*t41*t46-6336.0*t12*t24*t33*t4*t41)+3.0*t20*
        t21*t4*t71*grada+11.25*t19*t21*t4*t43*t52*grada-4.5*t19*t21*
        t28*t4*t49*grada-4.0*t11*t20*t4*t49*grada+3.0*t19*t20*t21*
        t48*grada+8.0*t10*t11*t48*grada+6.0*t11*t28*t4*t43*grada-9.0*
        t21*t25*t28*t4*t42*grada+6.0*t20*t21*t24*t42*grada-9.0*t19*
        t21*t24*t25*t28*grada-8.0*t11*t20*t24*t25*grada-6.5625*t1*
        t19*t21*t4*t72*t78-0.75*t1*t19*t21*t28*t4*t77-0.666666666666667*
        t1*t11*t20*t4*t77+0.5*t1*t19*t20*t21*t76+1.333333333333333*
        t1*t10*t11*t76-2.5*t1*t11*t4*t52*t72-2.25*t1*t21*t25*t28*t4*
        t71+1.5*t1*t20*t21*t24*t71+5.625*t1*t19*t21*t25*t4*t49*t52+
        5.625*t1*t21*t4*t42*t43*t52+5.625*t1*t19*t21*t24*t43*t52-2.25*
        t1*t21*t28*t4*t42*t49+3.0*t1*t11*t25*t28*t4*t49-2.25*t1*t19*
        t21*t24*t28*t49-2.0*t1*t11*t20*t24*t49+1.5*t1*t20*t21*t42*
        t48-2.25*t1*t19*t21*t25*t28*t48-2.0*t1*t11*t20*t25*t48+3.0*
        t1*t11*t24*t28*t43+3.0*t20*t21*t4*t42-4.5*t1*t21*t24*t25*t28*
        t42-4.5*t19*t21*t25*t28*t4-4.0*t11*t20*t25*t4+3.0*t19*t20*
        t21*t24+8.0*t10*t11*t24;
    res[13] = 0.5*t1*t20*t21*t4*(-4320.0*t5*t53*t7*t79/pow(rhoa,
        21.33333333333333)+4320.0*t5*t56*t81/pow(rhoa,18.66666666666667)+
        144.0*t4*t6*t76*t8*grada+432.0*t24*t48*t6*t8*grada+3456.0*
        t13*t16*t24*t4*t7*grada+18.0*t1*t4*t6*t8*t83+72.0*t1*t24*t6*
        t76*t8+54.0*t1*pow(t48,2.0)*t6*t8+216.0*t4*t48*t6*t8+216.0*
        t46*t6*t8+288.0*t13*t16*t22*t4*t7*t76-7200.0*t32*t5*t57*t74-
        3456.0*t24*t38*t4*t57*t74+3456.0*t24*t4*t54*t66*t7*t73+8640.0*
        t29*t5*t54*t7*t73-4824.0*t12*t30*t44*t5*t7+432.0*t13*t16*t5*
        t7-864.0*t30*t32*t4*t44*t48*t7+2160.0*t1*t13*t16*t4*t48*t7+
        864.0*t13*t16*t22*t24*t48*t7-864.0*t30*t32*t44*t46*t7+2160.0*
        t1*t13*t16*t46*t7-6336.0*t24*t30*t4*t40*t44*t7+2808.0*t1*t33*
        t45*t5+864.0*t12*t33*t4*t45*t48+864.0*t12*t33*t45*t46+5184.0*
        t22*t24*t33*t4*t45)+4.0*t20*t21*t4*t77*grada-8.0*t10*t21*t76*
        grada+15.0*t21*t4*t52*t72*grada-18.0*t21*t25*t28*t4*t49*grada+
        12.0*t20*t21*t24*t49*grada+12.0*t20*t21*t25*t48*grada-18.0*
        t21*t24*t28*t43*grada-1.0*t1*t10*t21*t83-6.5625*t1*t21*pow(t25,
        4.0)*t4*t78-3.0*t1*t21*t25*t28*t4*t77+2.0*t1*t20*t21*t24*t77+
        2.0*t1*t20*t21*t25*t76+7.5*t1*t21*t24*t52*t72+11.25*t1*t21*
        t4*t43*t49*t52-2.25*t1*t21*t28*t4*pow(t49,2.0)+3.0*t1*t20*
        t21*t48*t49+6.0*t20*t21*t4*t49-9.0*t1*t21*t24*t25*t28*t49-
        4.5*t1*t21*t28*t43*t48-12.0*t10*t21*t48-9.0*t21*t28*t4*t43+
        12.0*t20*t21*t24*t25;

}

static void
ft97bx_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[14];
 
    ft97bx_fourth_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];

    ds->df3000 += factor*res[5];
    ds->df2010 += factor*res[6];
    ds->df1020 += factor*res[7];
    ds->df0030 += factor*res[8];

    ds->df4000 += factor*res[9];
    ds->df3010 += factor*res[10];
    ds->df2020 += factor*res[11];
    ds->df1030 += factor*res[12];
    ds->df0040 += factor*res[13];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        ft97bx_fourth_helper(dp->rhob, dp->gradb, res);

    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

    ds->df0300 += factor*res[5];
    ds->df0201 += factor*res[6];
    ds->df0102 += factor*res[7];
    ds->df0003 += factor*res[8];

    ds->df0400 += factor*res[9];
    ds->df0301 += factor*res[10];
    ds->df0202 += factor*res[11];
    ds->df0103 += factor*res[12];
    ds->df0004 += factor*res[13];

}
